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Foreword 


This report was prepared by Dr. Reiner Rummel, Research Associate, 
Department of Geodetic Science, The Ohio State University (and now at the 
Deutsches Geodatisches Forschungsinstifut), Dr. Lars Sjoberg, Research 
Associate, Department of Geodetic Science, and Dr. Richard H. Rapp, Professor, 
Department of Geodetic Science. This work was supported under NASA Contract 
No. NAS 6-2484, The Ohio State University Research Foundation Project 
No. 7 839 04, This contract is administered through the NASA Wallops Flight 
Center, Wallops Island, Virginia 23337, with Mr. Ray Stanley, Project Scientist. ■ 

The three papers have been put together in a single report because of 
their common goal. Two of the papers are presented here for the first time, but 
one paper, by Rapp, was presented at the Spring American Geophysical Union 
Meeting in 1976.. 
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The Determination of Gravity Anomalies from Geo id Heights 
Using the Inverse Stokes' Formula 

by 

Reiner Rummel 



Abstract 


A numerical method for the determination of gravity anomalies from 
geoid heights is described using the inverse Stokes formula. This discrete form 
of the inverse Stokes formula applies a numerical integration over the azimuth 
and an integration over a cubic interpolatory spline function which approximates 
the step function obtained from the numerical integration. The main disadvantage 
of the procedure is the lack of a reliable error measure. 

The method was applied on geoid heights derived from GEOS-3 
altimeter measurements in the calibration area of the GEOS-3 satellite. The 
results for the estimation of several 5°, 2°, and 1° mean gravity anomalies 
are satisfactory when using a regularization cap with a spherical radius of 0. 1° 
and a reference field as e. g. obtained from the GEM 7 potential coefficients up 
to degree 16. 
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1. Introduction 


The anomalous potential T of the earth’s gravity field at a point P on 
the geoid is linked to the gravity anomaly at the same point via: 

<D . (-& + 7f) T( ' p ) 

which is called. the fundamental equation of physical geodesy, (c.f, Heiskanen and 
Moritz, 1967, p. 86) . In equation (1), y is the normal gravity and h the elevation 
along the normal. The zero order term of the potential is assumed to be zero. 

In spherical approximation and using Bruns formula to obtain the geoid height N 
instead of the anomalous potential equation (1) becomes: 

(2) Ag(P)= ~(^+ |-) GN(P) 


where R . . . mean radius of the earth, and 
G . . . mean gravity over the earth. 

This equation connects the gravity anomaly Ag with the geoid height N 
and has therefore to be the central equation for all further considerations. 

From the boundary condition (2) the solution of the classical problem, 
the determination of geoid heights from gravity anomalies, is derived the well- 
known Stokes' equation: 


(3) N (P) = ~~ J S Ag <Q) d <7 q 

a 

where S (<f>p q ) . .. Stokes' function for the spherical distance 
$pq . . . . . . between the points P and Q, and 
.cr unit sphere. 


For the solution of the inverse problem, i. e. the determination of the gravity 
anomaly field from the geoid height, we may use equation (3), too. It is then an . 
integral equation of the first kind. It will be an improperly posed problem in the 
sense that the unknown. gravity anomaly function is dis continuously dependent on 
the given function- N (P). One has to derive a regularized (here better: stabilized) 
solution of the integral equation (3). One way is to discretize equation (3) and to 
solve for mean values ,of a certain block size. Choosing this possibility even an 
overdetermined system of linear equations can be produced in order to find an 
approximate solution. For details, see (Gopalapillai, 1974), One may also aim 
for a global approximation to the gravity anomaly field by minimizing: 
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(4) min [||n (P), b s - S Ag (Q) ||^ + a 2 || Ag (Q) II* } 

Q > ^ n co co 

1 ■ CO 1 1 

The Stokes equation (3) is written in (4) in operator form. With expression (4) we 
minimize on one hand the discrepancy vector between the n given geoid heights 
and the corresponding geoid heights derived from the integration over the 
approximate gravity anomaly field Ag(Q), where ,T ~ " indicates "estimate". D 
is the error variance-covariance matrix with dimension (n x n). On the other 
hand, (4) tries to obtain ihe smoothest approximation Ag (Q), smoothest in the 
sense of minimum norm, where C is the global covariance function of the gravity 
anomaly field, e.g. the one derived by Tseheming and Rapp (1974). In expression 

(4) a is the regularization - (stabilization) - parameter that balances the first 
term against the second. The global character of the approximation should be 
indicated by the dimens ion 00 . The solution of the minimum problem (4) is; 

(5) Ag (Q) = C S T (S C S T + a D) -1 N (P) ob s 

co co co n co co n h 

1 co n co co n n 1 


which is for a - 1 identical to file least-squares collocation solution. It was 
successfully applied to this problem in (Rapp, 1974) and (Rummel and Rapp, 1977). 

Departing from the boundary condition (3) there exists alternatively to 
the above shown integral equation approach a direct inverse formulation of our . 
problem. It provides a direct means for the determination of gravity anomalies 
from geoid heights and was first given in (Molodenskii et al. , 1962):' 


( 6 ) 


Ag (P) = - N (P) 


_G_ P N(Q)-N(P) dap 
16 tt R J sin 3 jLP£ 

CT 2 


The derivation of equation (6) from the boundary condition (2) will be given in the 
next section. Equation (6) could also be considered as a solution to the Stokes 
formula (3) posed as an integral equation. Equation (6) again poses a regularization 
problem since for lim 4>pq 0 a small error d in N (Q) will cause a large error € 
in the determined gravity anomaly. 

This report will describe the discrete treatment of equation (6) given 
geoid heights as derived from altimeter observations. Results for estimated 
gravity anomalies computed from geoid heights will be presented as obtained 
from GEOS-3 altimeter data in the calibration area of GEOS-3. 
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Before we continue to discuss the numerical treatment of equation (6) 
it should be mentioned that gravity anomaly recovery using the integral equation 
approach, equation (3), as well as the inverse Stokes formula, equation (6), may 
also be tackled in the spectral domain. .This means, development into spherical 
harmonics for a global solution on the sphere. . The spectral form of equation (3) 
and of equation (6.) .are identical when expressed in spherical harmonics, namely* 


(7) g nD = ”• (n- 1) t„ B 

K 

where 


gnu'... spherical harmonic coefficient of the gravity anomal 
field of degree n and order m, 
t n0 ... spherical harmonic coefficient of the geoid height. 

The regularized form derived from equation (5) takes the form: 


(8) 




oba , n a 


with 


2 

a 8 variance of the observation error 

c n ....... i degree variance of the gravity anomaly- field 



square ratio of the radius of a chosen Bjerhammar 
sphere to the mean radius of the earth 


For the discrete form of equations (3) or (6), the analysis in the 
spectral domain means eigenvalue analysis of the corresponding system of 
linear equations. This will be more or less done for the discrete treatment • 
of equation (6) as shown in the sequel. Finally, when using a planar approxi- 
mation of equations (3) and (6), a Fourier analysis could be performed. 

Again, the spectral form of the planar of equations (3) and (6) will be identical 
in the spectral domain. 


2. Formulation of the Inverse.,Stokes Formula for Discrete -Data 

As already mentioned, we will shortly derive the inverse Stokes 
formula from the boundary value condition (2) which is in essence a spherical 
approximation of equation (1): 
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( 2 ) 


Ag(E) ~-G (®) + J. N(P) ) 


Provided N (Q) is twice differentiable for lim Q ~’P the radial derivative 5'N(P) 

B r ' 

becomes, using the Poisson integral (Heiskanen and Moritz, 1967, p. 38)‘: 


( 9 ) 


M B .i N(P) + r 

dr R K ’ 2rr J 


d Qq 

^o(P-Q) 3 


with i 0 . . . spatial distance between P and Q. Inserting equation (9) into 
equation (2), we obtain with: 


i&PQ 

& 0 (P- Q) = 2 R sin — - 

2 


the inverse Stokes formula 

( 6 ) 




_ sm — 
a 2 


The geoid heights derived from altimeter measurements are not given as a 
continuous function covering the sphere that approximates the earth as it would 
be required when using equation (6). Instead, there is only available: 

a finite amount of data 

arranged in arc segments with a certain along track and cross track spacing, and 

— disturbed by a certain amount of error. 

Thus we have to aim for a discrete solution of equation (6). The quality of the 
discrete formulation of the problem will entirely depend on the quality with which 
the radial derivative 9N (P) in equation (9) can be approximated from discrete data. 

& r 

The numerical treatment of the type of singular integral appearing on 
the right side of equation (9) is discussed in (Meissl, 1971). The procedure described 
there for the regularization of the integral cannot be applied in our case, not only 
because of the large computational effort that is necessary for its use, but mainly 
because of a lack of data dense enough to make this procedure applicable. 
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Integration Over the Azimuth 


Separating the terms dependent only of the spherical distance ip from 
those also dependent of the azimuth, equation (6) .is rewritten to: 

Ag (P) f (N (P) ♦ ^ J _/ <N(9>-.N(P» dor, d*,) 

tb ' 2 a * 

with 

sin ib = 2 sin cos 

2 2 

and with 

(10) • J (N(Q)-N(P» dtt Q = 2 tt'(N(^)-N(P)= 2tr 6 N(P, 0pq) 

a 

we obtain 

(11) Ag (P) = - -g (n(P) + | J 

The integration over the azimuth can be considered as an averaging- process that 
yields a function of' f>N (P, ipp q) for die computation point P only dependent of the 
spherical distance ip . As we know the accuracy, of such an-averaging process will 
be hardly affected by few bad data as long as many data enter the computation. 

For the numerical integration, the infinitesimal increment d ib with spherical 
radius if)q from the computation point is replaced by the finite increment Aip, 
with e. g. Aip = 0. 5°. These circular concentric zones are subdivided into four 
subblocks, as shown in Figure 1 with a certain arrangement of ground tracks 
of given geo id heights. A subdivision by four is chosen to obtain a representative 
mean value for each circular zone from the four subblocks even for a nonhomogeneous 
coverage of .the circular zone by altimeter tracks. A higher subdivision could not 
warrant that each subblock is covered by an altimeter track. 

First we average the given geoid heights in the four subblocks, then 
the four mean values are averaged to provide a mean value of the circular zone. 

This procedure yields a step function with step width Aip with the mean values of 
all circular zones which approximates the function 6N(P, $po) of equation (10), as 
shown in Figure 2. 


cos 


A 

2 


sin 


Z A 


3N(P, Ippq) d 0 q) 
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Figure 1: The templete of circular zones around the computation 
point P : As an example, all geoid heights (marked with ”o !l ) 
inside segment "a" of circular zone 3 will be averaged. 
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Uo 


I 



I 


Figure 2: Step function derived from the integration over the azimuths 
and cubic spline function approximating the step function. The 
integration over 0 is carried out from ipa to 0 na x . 
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Unfortunately, in the vicinity of the singularity point P, the 
circular zones are so small that only few geo id heights will be given inside them 
to be used for the average process. This may result in large errors for the 
mean values and may cause large errors in the determination of the gravity 
anomalies. A stabilization is therefore necessary to maintain sufficient accuracy 
at the expense of resolution, as will be described later. 


Integration Over the Spherical Distance 

We approximate the step function by an interpolatory polynomial 
function that attains the mean value of the i-th step (i.e. the mean value of the 
i-th circular zone with inner radius 0 t - A0 and outer radius 0 * + A0 ) at 

2 2 

the midpoint of the i-th step with spherical distance 0 i from the computation 
point p, as shown in Figure Two. This means the. polynomial function: 


(12) 


where 


P (0) = T aj (ip) 


j = i 


aj polynomial coefficient, and 

(0) . . . base functions 

fulfills the condition 


(13) 


6N (P,0i) = p(0i) 


We insert expression (12) into equation (11) and obtain: 


(14) 


Ag(P)=-Il(N(P)+ ^ J 

■tv 


J = 1 


0 


0 

cos 2 

sin 3 ! 

2 


(0) d 0 Q ) 


Since only the coefficients of the interpolatory series depend on the values of 
the step function the integration over 0 can be performed once for all with highest 
possible precision: 


(15) 


h 



0 


cos 2 



2 


(0) d 0 q 
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The gravity anomaly determination becomes, using equations (14) and (15): 


(16) 


Ag (P) = - (N<P) + aj Ij ) 

J = i 


For the interpolatory function a cubic interpolatory spline function is chosen.. 
Assuming n function values 6 n (P, 0 t ) to be given, it has the form; 


' (17) 




with 


+ ( 6N(P,» m ) M,^A <l> ) (# _ #0 + 


n - 1 


A0 =■ 0 i4 .i - 0i , and 

Mi = p" (0j), i. e. the second derivative of 

P(0) at 0 = 0i 


The advantage of using cubic splines is that the- coefficients Mi can be derived • 
from a simple recursion formula, (cf. Shampine and Allen, 1973, p. 54). For the 
unique definition of p (0) we define M 0 = M„.i = 0. The function p (0) defined that 
way fulfills the condition of being twice differentiable for lim 0-* 0, which is 
required for the function 6 N(P, 0) when using equation (9). 


define; 


Employing the basic principle of equation (16) and scaling A0= 1, we 


0i+i 


(18a) 


if 1 - cc 

= jg J- •t(0 1 + l-0) 3 - (01+1-0)] - 


0 

cos 2 


' 0r 


sin 2 l 
2 


d 0 


0 


< 18b > ^ - h I 


i+i 


0: 


0 

„ COS 7T 

U01-0.) 3 - (0.1-0)] — r- d 0 

sin JO. 

2 
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<3$ 








0? 



„0 


1+1 


•j p 1 COS o 

Is.i = xJ — XT d^,and ( 

0 t sir. ^ 


0 


(18c) 


sin r 
2 


■iv.yri .■<; : 


0 


(18d) 


U, l 


01+1 COS-r- 

f j <**- *) r4- d 0 

4 J a ?» 


03 


sin 2 ! 
2 


The integration is carried out numerically in extremely small steps by Gaussian 
quadrature. Inserting equations (18, a - d) into equation (11), the estimation formula 
for a gravity anomaly at a point P becomes: 

_n— 1 

(19) Ag (P) “ “ ir (n(P) + / [Mi Ii f i - Mi+i la, i + 

i=o 


+ 6n(P,0i) I 3 ,i - 6N(P, 0 1+ i) 14, i]) 

In equation (11) the integration ranges from 0=0° to 0= 180° which would require 
data given over the entire globe. Actually, because of the rapid decrease of the 
integral kernel 1 in the inverse Stokes formula,, accepting a small error, 

sin 3 ± 

2 

the integration can be truncated at a suitable 0 aaX . Gopalapillai (1974) investigated 
the error of trun cation . From his results we conclude that an integration up to 
0 max = 10° suffices at the present accuracy level for the geoid heights. - This means 
that we accept an error of a. tr ±0.3 mgal for the estimated gravity anomalies when 
a reference field up to. degree n,,,^ 16 is introduced and an error of CT tr = ± 3 mgal 
for a reference field with n asx = 4. 

The integration from 0= 0° causes, as already noted, stability problems. 
At the critical point 0= 0°, the integral kernel goes to infinity with the strength 
1 

sin 3 j£ » whereas for the difference 6 n (P, 0 p q) we know only it has to be twice 
2 

differentiable. A small error 6 in the geoid height difference close- to the 
singularity point will therefore blow up to a large error e in the estimated gravity 
anomaly. Because of the lack of dense data, the regularization method proposed 
by Meissl (1971) cannot be employed in our case, so we introduce a very pragmatic 
method of regularization. We postulate that the geoid height difference 6N(P, 0pq) 
is zero in a small spherical cap around P with spherical radius 0 0 . In our 
computations, the radius of the "regularization cap" was varied between 0o = 0.001° 
and , 00 = 0 . 1 °, which corresponds to a distance of approximately 0.2 km to 22 km. 
The. integration is thus actually carried out from 0= 0 O to 0 nax . 
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The inverse Stokes formula, equation (6), provides point values at 
Pea. What we are interested in are mean values- of e.'g. block size l°x 1° or 
2°x 2°. The correct determination of mean values of a block would require an 
integration over the block or. in other words, the geoid height function given 
continuously inside this block. The procedure used in our computations first 
estimates all point values of gravity anomalies inside a block for which discrete 
geoid heights of an arc segment are given, using equation (19). Then these 
point estimates are simply averaged to give the mean gravity anomaly of the block. 

As a severe .drawback of the use of the inverse Stokes formula in discrete 
form, we have to accept the fact that there is no easy way to find a reliable error 
estimate to the estimated mean gravity anomalies. 


3. Numerical Results 

The data used in the computations were geoid heights as derived from 
GEOS-3 altimeter measurements in the calibration area of the GEOS-3 satellite. 
The distribution of the source data and the method of the geoid height determination 
are described in (Hummel and Rapp, 1977). We estimated mean gravity anomalies 
for the same blocks where the mean gravity anomaly recovery was performed using 
the least squares collocation technique, equation (5), and also described in (Ibid. , 
p. 82). These are three 5° equal area anomalies, two 2°x 2° anomalies, and two 
l°x 1° anomalies, as shown in Figure Three. The location of these blocks is 
chosen to maximize the number of altimeter tracks given at that time in and around 
the blocks. 


The computations were carried out with the geoid heights referred to 
an ellipsoidal reference field and with the geoid heights referred to a reference 
field up to degree n max = 16 computed from the coefficients of the Goddard Earth 
Model 7, (Wagner et al. , 1975), using: 

n _ _ 

(20) N(P)(SM7=~~y ) (C^ cos mA p+ Sno sin m A p) P nn (sin ( 0 p) 

7R Li L> 

n=2 m=0 


for the computation of the reference geoid heights. All terms of equation (20) are 
explained in (Rapp and Rummel, 1975, p. 2). We have then to add to the estimated 
residual gravity anomaly a reference gravity anomaly also computed from the same 
set of potential coefficients up to degree n nax using equations (20) and (7). 

The results of the gravity anomaly recovery using the discrete form of 
the inverse Stokes formula are shown in Table One, Columns 7 to 11. They are 
given for different regularization caps ($ 0 -= 0.001°, 0.01°, 0. 1°) and for an 
ellipsoidal as well as for a low degree (n nax = 16) reference field. In column 12 
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Figure 3. Location of Test Blocks and Related Altimeter Tracks 
in the Anomaly Estimation Area 
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the number of point values computed to obtain the mean gravity anomalies are 
given. Table One gives also the results of the mean gravity anomaly recovery 
when using the least squares collocation method, taken from (Rummel and 
Rapp,' 1977, Table Five). In column 4 the known terrestrial mean gravity 
anomalies are shown for the same blocks; in column 13 the mean gravity 
anomalies as derived solely from the potential coefficients of the GEM 7 
earth model up to n BaX = 16 are given. 


4. Conclusions 

The results of Table One for the estimation of mean gravity anomalies 
using the inverse Stokes formula show: 

that a regularization cap witlra too small spherical radius ip 0 (ip 0 ~ 0. 001° and 

also ip o = 0,01°X yields too unstable .results. Only for a ip Q = 0,.1° stable estimates 
could be obtained. 

only the- use of the reference field (here n- ttax = 16) could provide satisfactory 

results, whereas the gravity anomaly estimates using the ellipsoidal -reference 
field were to far off in comparison with the terrestrial- gravity anomalies. 

Using the reference field (GEM 7, n Bax = 16) and a regularization cap with )p o= 0.1° 
the presented method for the .use of a discrete form of the inverse Stokes formula 
produced satisfactory results. The results, as shown- in column 11 of Table One, 
are in good agreement with the given terrestrial gravity anomalies (column 4) and 
with the results obtained from collocation (columns 5 and 6). Thus, the described 
method can be considered a valuable alternative to the least squares collocation 
method. 

Its advantage in comparison with collocation will.be that it needs less 
computer time; no system of linear equations has to be solved. The - d rawbacks are: 

there is ho easy way' to* find a reliable error estimate to the determined mean 

gravity anomalies'. 

the method of regularization is from -the theoretical point of view hardly 

satisfactory. 

- ■ I 

For a further improvement -of the presented method^ one should therefore investigate 
more appropriate buU'still economical' methods of regularization. In addition, 
instead of the type of spline function used for these computations, one should use 
a mean value reproducing spline, as described in (Siinkel, 1975) which, when 
integrated over the intervals A ip , yields exactly the ordinate value of the step 
function. More stable results could be obtained using an adjusting spline function 
with less nodes than given steps of the step function. Finally, we would expect to 
obtain better gravity anomaly recovery results given denser geoid height data without 
significantly increasing the computer time. - ■ 
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Table One: Results of the estimation of mean gravity anomalies using the inverse Stokes formula 


O 

I 


O. 


1 

2 • 

3 

4 

5 

6 

■a 

8 

9 

10 

11 

12 

13 

block 

size 

coordinates 

terr. 

collo- 

collo- 

inv. 

inv. 

inv. 

inv. 

inv. 


Ag 




Ag 

cation* 

cation* 

Stokes 

Stokes 

Stokes 

Stokes 

Stokes 


from 



<P 


(elllp.) 

(GEM 7 

(ell ip.) 

(GEM 7 

(ellip.) 

(GEM 7 

(GEM 7 


(GEM 7 






(12,12)) 


(16, 16)) 

>l>0~ 

(16,16)) 

(16,16)) 

| 

(16,16)) 



A 




0.001° 


ESESfl 

^0- 

>l> 0= 

■ 









0.001° 


0.01° 

0,1° 

■ 




45°- 40° 











D 

5° E. A, 

299°- 30 6° 

-17 ± 3 

-16 ± 6 

-20 ±3 

-23 

• -31 

-13 

-26 

-21 

113 

-15 



40°- 35° 











B 

5° E. A. 

291°-297° 

-24 ± 2 

-24 ± 5 

-25 ±3 

-10 

-31 

- 7 

-27 

-24 

93 

-16 



35°— 30° 











c 

5° E. A. 

295°- 301° 

- 8 ± 3 

- 3 ± 6 

-11 ±3 

12 

14 

11 

-14 

-15 

97 

-23 



40°- 38° 











D 

2° x 2° 

291°-293° 

-31 ± 7 

-33 ± 7 

-32 ±5 

- 9 

-41 

-10 

-39 

-36 

26 

-22 



34°- 32° 











E 

2° x 2° 


- 8 -fc 7 

- 2± 7 

- 8 ±5 

11 

- 6 

10 

- 7 

- 5 

23 

-13 


- 

296°-298° 











/ 


39°- 38° 











F 

- 1° x 1° 

291°- 29 2° 

-25 ± 14 

-30 ±11 

-29 ±9 

-27 

-38 

-21 

-33 

-29 

8 

-12 



34°- 33° 











G 

1° x 1° 

296°-297° 

- 8 ±15 

- 3 ±10 

- 9 ±8 

6 

-13 

7 

-12 

-14 

9 

' -22 


* from (Rummel and Rapp, 1977) 
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Abstract 


A method is outlined for the prediction of mean free -air gravity 
anomalies from altimeter data along one arc at a time. The method may be 
regarded as least squares collocation performed in the frequency domain. In 
this way a considerable gain in the necessary computer time is achieved. The 
final prediction of a block mean anomaly is estimated by the mean value of all 
separate estimates from different arcs. A simple formula is derived for 
estimating the prediction errors . . 

The method is tested for the prediction of 5° equal area anomalies 
and 1° x 1° anomalies in the GEOS-3 Calibration Area and in the Philippines 
Area. For example, approximately 12800 observations in 85 arc segments were 
processed for the prediction of 6 5° anomalies in the Calibration Area. The 
total computer time was 2. 1 minutes and the RMS difference between predicted 
and terrestrial anomalies was 4 mgal. In the same area a subset of 21 1° x 1° 
anomalies was predicted from 41 are segments (approximately 6500 observations) 
with a corresponding RMS difference of 5 mgal. 

Due to the simplicity of the method, all available altimeter arcs in 
an area can easily be included. Additional information may in a simple way be 
taken into account in a posterior prediction. 
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1. Introduction 


The method of least squares filtering of data measured at discrete 
points along a profile was presented by H. Moritz (1967). In 1969, Moritz 
extended the method to the prediction in the area between observed profiles. 

The latter technique suffers from the disadvantage of requiring parallel 
profiles of observations. However, Moritz* idea of predicting in points not 
on a profile may be generalized to the prediction of gravity anomalies from 
altimeter data given along one profile (arc) at a time. Then the final estimates 
are determined by the mean value estimates provided by all the arcs in the area. 

The method may be regarded as least squares collocation performed 
in the frequency domain (instead of in the covariance domain). Subsequently, 
the covariance functions are transformed to the corresponding spectra. The 
autocovariance matrix is transformed to diagonal form, which considerably 
reduces the time of computation. Thus the very large amount of available 
altimeter information can be handled in a simple way. 


2. Theoretical Background 

Let us assume that a stochastic process x is measured at n points 
( discrete case )'. The process y related to x can then be estimated linearly by 
the observations x t in the following way; 

n 

(2. 1) yp = ^ hpi x t 

i=i 


where h PI are the weights of the observations. The prediction error is given by: 


n 



1=1 


The linear least squares prediction problem is to estimate' the weights in such a 
way that the RMS prediction error is a minimum. We assume that the statistical 
expectations of x and y are 0 : 


( 2 . 2 ) 


E { x ] = 0 and E { y } = 0 
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For each set of weights the prediction variance becomes: 


(2.3) 


n 

mp- = E {ef } = C(yp,y.' P ) -2 Y h Pl C (yp,x t ) +■ 

1= i 


n a 

+ 1! I! k pi hpj c ( Xi * Xj ^ 
1=14=1 


where 


C (y P ,yp) = E f yf } 

C (yp.xj) = E {y? xj} 
C (xj, xj) = E {Xi xj 


variance of y p 

cross-covariance of yp and xi 
auto-covariance of observations 


The minimum of mp is obtained by differentiation of (2. 3) with respect to each 
of the weights h P * , which gives the (generalized) Wiener-Hopf equations in the 
discrete case: 


(2.4) £ h Pi C(Xi,x d ) = C(y P , xj) ; j = 1, 2, ..., n 

1=1 

or with obvious matrix notations 

(2.4) h C** = c yx 
with the solution 

( 2 . 5) h = c y x C xx 


By inserting the solution for h into (2. 3), the least squares prediction variance 
becomes (with matrix notations) : 


( 2 . 6 ) 


m p = C(y P ,y P ) - c yx Cxx 1 c yx T 
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Let us now study the corresponding continuous problem. Suppose 
that the stochastic process x is recorded continuously along a profile, and that 
the process is stationary (such an assumption is normally imposed also in the 
practical application with discrete data as outlined above). Formula (2.1) is 
then defined; 


(2.7) 


ym - J 


CO 

h (t - tj x(t) dt'= J h (T) x (t-T) d t 

— GO 


and the prediction variance (2.3) becomes 


( 2 . 8 ) 


m = cr v 


w 00 CO 

- 2 J h (a) Cyx (a) d a + J J h(a) h(j8) C xx (a - 0)da d/5 


where of is the variance of y. In this case the least squares prediction problem 
is solved by the following Wiener-Hopf integral equation (Papoulis, 1965, p. 345); 


(2.9) 


CO 

• J h(tt) Cxx (T- a) da = Cyx ( T ) 

— 03 


The weight function h (Of) is most conveniently obtained from the Fourier transform 
of (2. 9) . Let us use the notations : 


H(CO) = K(T) , S xx (OJ) = C XX (T) , S, yx (CO) = C yx (T) 


where X denotes the Fourier transform of X, i.e. : 

CO 

(2.10) X(to) = J X(T)e' iUT dT 

— CO 


Using the convolution theorem (Papoulis, ibid. , p. 159) the Fourier transform of 
(2.9) becomes: 


H(£ 0 ) S xx (co) = S yx (CO) 
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Thus 


( 2 . 11 ) 


H (to) = Sy X (CO) / S** (to) 


Finally, the weight function h(T) is given by the inverse Fourier transform of H(tO): 


( 2 . 12 ) 


00 

1 P itOT „ 

h < T ) = 2 ^ J H (to) e dto 

— CO 


The prediction variance along the profile may be determined directly from the 
spectrum of (2, 8), which becomes; 


CO 

(2.13) m 2 = — J [Syy(to)- H(w) Syx (to)] d to 

— CO 

The basic formulae for this study are (2.7), (2.11), (2.12) and (2.13). In the 
numerical applications, the integrals are approximated by summations and the 
integration limits are truncated at finite values. 


3. Prediction of Anomalies Along Arcs 

In this section we are going to present the formulae to be used in a 
practical application of the Fourier transforms for the estimation of free-air 
point gravity anomalies along the arc from an (adjusted) altimeter data profile. 

The corresponding formulae for least squares filtering are given in Moritz 
(1967). By using a low degree spherical harmonic expansion as the reference 
field, the remaining altimeter residuals can be regarded as regional and local 
variations with mean value 0 { cf. (2.2)} , The purpose of the method is to 
improve the gravity anomalies provided by the spherical harmonic expans ion 
by adding the information obtained from altimeter data. Rigorously, the method 
requires that the data have equal spacing. This will be used as an approximation 
for each arc . 

The point gravity anomaly Ag k is predicted from the residual altimeter 
height (geoidal undulation) X = N - N by {cf. (2.1)} : 

(3.1) A& k = Ag k + At £ h j • (Xk-j + Xk +J ) 

o 
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where 


N = reference geo id undulation (in meters) 

Ag n = reference gravity anomaly (in mgals) 

At = spacing of data (in meters) 

m = selected number of discrete intervals of spectra 


At is determined from the spherical distance of the total arc segment and the, 
number (n) of observations : 


(3.2) 


At = r n 0/(n - 1) 


where r B is the mean earth radius in the area. Furthermore, the fundamental 
wavelength: 


(3.3) 


A co = n/m At 


is computed. The choice of m has to be based on the required resolution of the 
solution. It remains to determine the weights h j , which is done in the following 
way. 


Covariance functions. In Moritz (ibid. ) it is suggested to compute the covariance 
functions (Cnn> C^ s ,n and C^ s ^ 6 ) from observations. However, as the only kind of 
observations available in the present case is altimeter data, we will use the empirical 
covariance functions given by the subroutine COVA of Tscherning and Rapp (1974). 

The low degree harmonics included in the reference field are deleted in the 
covariance functions. Computer time is gained by storing the covariance functions 
in a discrete table and then determining the necessary values by interpolation 
between table values. The covariance functions are computed for: 


Cnn (4 At), C^ gN (A At) and C& g ^ (4 At); 4=0, 1 m 


The observation errors are taken into account in the autocovariance matrix Cnn , 
which is replaced by: 

-Cnn (0) + CTn 2 if 4 = 0 
(3.4) Cnn (4At) = ■{ 

l Cnn (4 At) if 4^0 
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where cr 2 is the mean square error of the observations along the arc: 



CTj being the standard error of the i- th observation. 

Spectra. The following Fourier transformations are performed; 

Cnn Snn , C A6 n S^g n and C A g A g “* S AK A g 

The numerical integration is as follows [cf. formula (2. 10)] : 

Raw values; 


B- 1 


(3.5) 


where 


S/ = At [C o + 2^ Cq cos (qnr/m) + C n cos (r tt) 


q= 1 


r = 0, 1, 2, ...» m 


Smoothed values: 


(3.6) 


So = “ (S o + Si) 

Sq = ~ (Sq -1 + 2 S q + Sq +1 ) ; q = 1, 2 m - 1 

S m = — (S B _i + S n ) 


System function . 


(3.7) 


Hr = S A g N (rAco)/ S NN (r Aco) ; r = 0,1, 2, . . . , m 
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Weight function. Inverse Fourier transform [See (2. 12)]: 


H (O)) - h (T) 


Raw values: 


(3.8) 

where 


a- 1 

hjj = [ Ho + 2 Hq cos (qirr/m) + H„ cos (Itt)] 

< 1=1 


& 0 j 2j » • • f 21 


Smoothed values: 


(3.9) 


ho - ^ (ho + hi) 

hq = ~ (hq-1 + hq' + hq+ 1) ; q = 1, 2, . . . , m - 1 

h„ — “ (h E _i + 2 h 0 ) 
h B +x ^ h E 


In. this way the necessary weights of (3. 1) are determined. 


Remark. In Moritz (Ibid., p. 49) the factor 1/4 of h 0 is incorrectly set to 1/2. 
For further details see the Appendix. 


4. Prediction of Mean Anomalies 


In Moritz (1969) a method is given for a least squares interpolation of 
data between parallel profiles. As the altimeter arcs are usually not parallel, this 
method cannot be used. However, the technique of predicting in points that are not 
on the profile of observations may be used for one arc at a time (extrapolation). 

This method requires that the shortest distance (a) from the profile to the point 
of computation is computed. Then the cross-covariance matrix used in the previous 
section is simply modified to: 
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(4.1) 


Ca?n ( t ) = Gas N (a) , a = /a 2 +^ 1 

where G is the output from subroutine COVA. The auto-covariance matrices 
Cnn and C^ Ag are not modified. Thus the auto-spectra Snn and S AgAg can be 
computed once and for all for each profile (arc). The cross-spectrum S AgN and the 
weight function h(T) must be recalculated for each new distance a. However, for 
a constant a the same weight function is used for repeated predictions. This fact 
allows for the computation of many points inside each block without significant increase 
of the cost. A further s im plification of the mean block predictions is obtained by 
calculating C AgN (r) as the average of three different a -values for each block: 

i 

(4.2) Cas M (t) = V G (^ j2 + (a+ ]yf)/ 3 

J=-i 

where a is the shortest distance between the block center and the arc, and y is 
approximately half the side of the block, y is selected a priori by the user. If 
the profile was parallel to one side of the block, y could be selected half of the 
perpendicular side without introducing any error. However, in practice, the 
profiles are generally not parallel to any side of the block and for any choice of 
y an approximation error is present. 



Figure 4, 1 . The cross-co variance function is applied as the average 
of the covariance functions among the point p on the profile and the points Qi - Q 3 
on Hie normal to the profile, y is approximately half the side of the block. The 
predictions are performed for all points P with P x ^ P ^ P2 , and the block mean ' 
prediction is represented by the mean value of all individual predictions' (in the 
shaded area). 
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5. Accuracy Estimation 


The variance of prediction of y along a profile from the observations 
x was given in formula (2. 12). In the practical application, this formula takes 
the form: 


(5. 1) 


m 2 = 


Ato 

2tt 


a — 1 

^E o + 2^ Eq + E n ^ 

(i = i 


where A to was defined in (3.3) and 


Eq = Syy (qAco) - Eq s yx (q Aw) ; q = 0, 1, 2, . . . , m 


Now we proceed to estimate the errors of mean gravity anomaly predictions. The 
error covariance (after prediction) between the points p< and P , are given by 
(Moritz, 1969. p. 41) ; 


(5.2) 


00 

i r icoT , 

CTlj (T) = 2^ J El} (C0) 6 d W 

~ CO 


where 


ElJ (<D) = S yy (tO) - H (tO) Syx (to) 
t = distance between P i and P j 


We define the mean value of the prediction y (t) over the interval t x £ t £ t a by: 

' ^3 

(5. 3) y = _L J y (t) d S 

^1 


where 


b = (t 3 - tx)/2 
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The variance of' y is then given by: 

oo 

(5.4) m 2 =^=J -E,j(W) F(CO)dCO 

— CO " s 


where 


F (co) = 


r 


l 


1 if co = 0 

[sin(b co.) /1b co 3 ] if 


co^ 0 


Proof: Let us denote the error of y(t) by € (t). The error of the. mean (y) is 
then given by: 

ta 

k J E < t)dt 

t.i 


and 


m 2 = E [€ 3 3= -±- 

4 b 2 


t s ta 

e{JJ €(t) e(t') dt dt' j- = 

tx ti 


1 

4t? 


ta ta 

J J CT ?J (t - 1) dt d.t' 

tl tx 


Inserting (5. 2) we obtain: 


<» 

m 2 = “ J E u (co) F (W) d co 


where 


F( “ ) = ip J I e '“ (t ~‘" )dt dt '= { 


t x 


T for co = 0 
(sin(bco)/bco) 2 for co^O 
Q. E. I). 
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Using the notation 


Jr = E i j (to r ) F (tor) 


We finally arrive at the following practical formula: 
(5.5) 

1 = 1 


m 


A to f 

2 tt \ 


-i 


J o + 2^T J q + J 


") 


As the cross-covariances that are used in the prediction of mean gravity anomalies 
are assumed to be mean values across each block perpendicular to the altimeter 
profile (cf. section 4), formula (5.5) may be applied directly for estimating the 
prediction error of the block mean value. In this case, the constant 2b is the 
side of the block. 

Special cases* 

a) b -» 0 implies F (to) -» 1 for all to and therefore in 3 -* m 2 , i. e. we arrive at the 

ease of point estimation. • - 

b) b -* 03 implies F (to) -» §(to) and m 3 -» Ey (0)/2n. ' 


6, Computations 

In all computations equidistant spacing of the data along each arc is 
assumed, and the mean earth radius (r s ) is set to 6370 km. First, the Fourier 
transform method described in section 3 was tested for least squares filtering of the 
data along some arcs. The filter is given by formula (3. 1) with Ag k and Ag k 
substituted by Nk and Nk and S^ sN of (3. 7) substituted by Snn • The GEM 7 spherical 
harmonic expansion to degree 16 was used as a reference. An example of this test 
is shown in Figure 6.1. A reasonable smoothing of the raw altimeter data is 
obtained. 


Second, mean gravity anomalies were predicted from altimeter data. 

The spherical harmonic coefficients -.of GEM 9 to degree 20 were used, in the 
reference field. The covariance output from subroutine COVA was stored in 
tables at 0.°1 intervals. Each profile was processed separately. The final mean 
anomaly of a block (Ag) was estimated by the weighted average of the outcomes 1 f * 
(Ag i) from several arcs in the area: 
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Figure 6.1. Filtering of 300 alti 
Fourier transforms 
observations is 0 . 6. 
spectra* 75. 



( 6 . 1 ) 


Ag = m 2 



m 


2 

i 


i = I 


where m 2 is the estimated variance of the mean: 


( 6 . 2 ) 




m 


-2 


i= i 


and m f are the estimated variances of predictions Ag t . 

Requirements of used arcs . Each segment of an arc with less than 100 sequential 
observations [with gaps less than 0?5 (0°3) for 5° (l°x l e ) anomalies] were rejected. 
The maximum accepted distance from an arc to the center of a block was 3° (if not 
specially mentioned) . For details of the adjusted altimeter data, we refer to Rapp 
(1977 b). 


The cross -covariance matrix C^gN was computed as the average of 3 
points perpendicular to the profile [see formula (4, 2)], The 5° (l°x 1°) mean 
anomalies were estimated as the averages of 21 (5) such sequential "mean-value 
points" predictions inside each block parallel to the profile. The area thus 
covered with point predictions is illustrated as the shaded area of Figure 4.1. 

Again we emphasize that for any choice of y and number of "mean-value points" 
inside a block for forming the block average, an approximation error occurs due 
to the inclusion of points outside the block and (or) the exclusion of points inside 
the block (Figure 4,1). However, this type of error is more or less present in 
any numerical integration and as long as the error does not become critical, the 
quality of the final prediction should be compared with the necessary computer 
time. In this respect, the selected method is favorable. 

The number of discrete intervals (m) in the numerical integrations of 
the spectra was set to 60 and 75 for 5° and l°x 1° anomaly predictions, respectively. 
The prediction results are compared with the terrestrial mean anomalies described 
in Rapp (1977a). 


6. 1 The Calibration Area 


Data : Adjusted altimeter data in the area 10° < <p < 50° and 275° < A < 310° (see 
Rapp, 1977b). 

In this area 6 5° equal area anomalies and 30 l°x 1° anomalies were 
estimated (all the 1° x 1° blocks are inside 5° block No. 403; see Rapp, 1977a). The 
predic tion results of the 5° blocks are shown in Tables 6.1 and 6.2. These compu- 
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tations include approximately 12800 observations in 85 arc segments, and the 
total computer time (IBM 370/168) was 2. 1 minutes. The RMS differences 
prediction- terrestrial and GEM 9 - terrestrial anomaly are 3. 7 and 7. 1 mgal, 
respectively. 

The prediction results of the 30 l°x 1° anomalies are given .in 
Tables 6. 3 - 6. 4. Approximately 6550 observations in 41 arc segments were 
processed, and the total computer time was 8 minutes. The RMS differences 
as above are 11. 5 and 13.8 mgal. However, as the RMS standard error of the 
terrestrial data is as large as 13.6 mgal, these differences might well be due 
to the uncertainty of the terrestrial gravity information alone . Subsequently, 
by excluding the 9 last blocks (with poor terrestrial information) from the 
comparison, the remaining subset of 21 blocks gives the following RMS differences: 


| Ag A - Agt | = 4.8 mgal and | Ag G EM - Agr 1 = 6.7 mgal 


Table 6. 1 


5° Equal Area Predictions in the Calibration Area 
Block Numbers and Terrestrial Anomalies (Agi ± mj) Refer to Rapp (1977a). 

Units? mgal 


Block 

Number 

Coordinates 

Predic 

. Means 

* 

KlA 

Spread 

No. 

of 

Arcs 

(n) 

Agr± m T 

GEM 

9 

Wtd. 

Direct 

342 

40°-35°, 284°-291° 

-23.7 

-24.7 

1.2 

3.8 

41 

-21.9 ±1.7 

-17.0 

343 

40°-35°, 291°-297° 

-20.0 

-20.1 

1.0 

3.5 

37 

-24.3 ±2.0 

-17.5 

402 

35°-30°, 283°-289° 

-28.0 

-29.9 

0.9 

4.7 

44 

-34.1 ±2.4 

-20.-6 

403 

35°-30°, 289°-295° 

-24.5 

-23.1 

1.0 

5.9 

44 

-26.8 ±2.5 

-27.7 

465 

30°-25°, 281°-287° 

-18.2 

-17.3 

0.9 

2.9 

41 

-15. 4 ±2. 8 

-19.2 

466 

30°-25°, 287°-293° 

-33.4 

-31.2 

1.0 

4.1 

34 

-29.8 ±2.8 

-35.9 


* estimated prediction error from altimeter data 
+ = [E (Agi - Ag) 2 / n]a , where Ag t = prediction from arc i. 
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Table 6.2 


Differences Between Predicted, Terrestrial and GEM 9 Anomalies 

in the Calibration Area 
Units : mgal 


Block 

No. 

Weighted 

Terr. 

Direct 

Terr. 

. GEM 9 
Terr. 

Weighted 
GEM 9 

342 

-1.8 

-2.8 

4.9 

-6.7 

343 

4.3 

4.2 

6.8 

-2.5 

402 

6.1 

4.2 

13.5 

-7.4 

403 

2.3 

3.7 

-0.9 

3.2 

465 

-2.8 

-1.9 

-3.8 

1.0 

466 

-3.6 

-1.4 

-6.1 

2.5 

RMS: 

3.8 

3.2 

7.1 

4.5 
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Table 6. 3 


30 1° x 1° Mean Anomaly Predictions from Altimeter Data 
in the GEOS-3 Calibration Area. 

Agr and mr According to Rapp (1977a) 

Units: mgal 


<p°, A 0 
(NW comer) 

Predie. 

Means 

m A * 

Spread 

#of 

Arcs 

AgT 

± mi 

GEM 

9 

Wtd. 

Direct 

35 

289 

-28.0 

-28.0 • 

3.1 

6.3 

19 

' -33 

6 

-22.9 

35 

290 

-27.2 

-26.3 

3.0 

5.3 

20 

-24 

7 ‘ 

* -23.6 

35 

291 

-25. 2 

-24.4 

3.0 

4..6 

19 

-29 

2 

-23.8 

35 

292 

-22.7 

-21.1 

2.8 

5.1 

25 

-30 

2 

-23.6 

35 

293 

-19.3 

-17.3 

. 2.7 

5.8 

29 

-24 

3 

-23.1 

35 

294 

-14.4 

-13.5 

3.2 

6.4 

26 

-17 

11 

-22.2 

34 

289 

-29. 6 

-30.6 

3.5 

7.6 

18 

-32 

14 

-25.1 

34 

290 

-26.1 

-26.9 

3.4 

6.6 

19 

-31 

11 

-25.8 

34 

291 

-24.1 

-23.4 

2.9 

5.9 

24 

-26 

15 

-26.1 

34 

292 

-21.7 

-19.7 

2.8 

6.6 

27 

-25 

10 

-25.8 

34 

293 

-17.9 

-14.8 

2.9 

9.3 

27 

-25 

11 

-25.2 

34 

294 

-11.9 

- 9.6 

3.2 

10.9 

25 

-16 

15 

-24.2 

33 

289 

-32. 8 

-32.9 

3.4 

7.1 

20 

-38 

16 

-27.4 

33 

290 

-28. 6 

-28.5 

3.2 

6.5 

22 

-32 

18 

-28.2 

33 

291 

-24. 6 

-24.2 

3.2 

6.7 

24 

-28 

12 

-28.4 

33 

292 

-20. 6 

-19.7 

3.2 

7.6 

26 

-30 

14 

-28.2 

33 

293 

-15. 5 

-13.4 

3.2 

10.4 

24 

-25 

14 

-27.4 

33 

294 

- 9.4 

- 5.6 

3.2 

19.8 

25 

- 9 

6 

-26.3 

32 

289 

-36. 8 

-36.0 

3.3 

6.7 

19 

-40 

19 

-29.6 

32 

290 

-32.3 

-31.0 

3.2 

7.1 

24 

-35 

19 

-30.5 

32 

291 

-27.2 

-26.2 

3.1 

7.5 

27 

-26 

13 

-30.8 

32 

292 

-22.3 

-21.2 

3.2 

8.5 

26 

-41 

22 

-30.5 

32 

293 

-18.8 

-17.3 

3.3 

8.8 

23 

-23 

17 

-29.8 

32 

294 

-12.2 

-12.4 

3.7 

8.7 

21 

-15 

14 

-28.5 

31 

289 

-36.1 

-35.2 

3.0 

7.0 

23 

-11 

15 

-31.8 

31 

290 

-33.3 

-31.6 

3.2 

7.6 

26 

4 

15 

-32. 7 

31 

291 

-29.1 

-26.3 

3.3 

8.8 

25 

- 2 

15 

-33.1 

31 

292 

-23.3 

-21.5 

3.2 

10.5 

25 

- 6 

15 

-32.9 

31 

293 

-17.1 

-16.3 

3.6 

9.8 

21 

-10 

15 

-32.0 

31 

294 

-13.7 

-14.6 

3.9 

8.4 

19 

- 8 

14 

-30.7 


* * estimated prediction error from altimeter data 

+ = [E(Ag i - Ag) 2 / n]£ , where Ag t = prediction from arc i 
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Table 6. 4 


Differences Between Altimeter (Weighted and Direct Means), 
Terrestrial and GEM 9 Anomalies According to Table 6. 3 

Units; mgal 


(NW corner) 

W - Terr. 

D- Terr. 

GEM 

Terr. 

W- GEM 

35, 

289 

5.0 

5.0 

10.1 

-5.1 

35, 

290 

- 3.2 

- 2.3 

0.4 

-3.6 

35, 

291 

3.8 

4.6 

5.2 

-1.4 

35, 

292 

7.3 

8.9 

6.4 

0.9 

35, 

293 

4.7 

6.7 

0.9 

3.8 

35, 

294 

2.6 

3.5 

- 5.2 

7.8 

34, 

289 

2.4 

1.4 

6.9 

-4.5 

34, 


4.9 

4.1 

5.2 

-0.3 

34, 

291 

1.9 

2.6 

- 0.1 

2.0 

34, 

292 

3.3 

5.3 

- 0.8 

4.1 

34, 

293 

7.1 

10.2 

- 0.2 

7.3 

34, 

294 

4.1 

6.4 

- 8.2 

12.3 

33, 

289 

5.2 

5.1 

10.6 

-5.4 

33, 

290 

3.4 

3.5 

3.8 

-0.4 

33, 

291 

3.4 

3.8 

- 0,4 

3.8 

33, 

292 

9.4 

10.3 

1.8 

7.6 

33, 

293 

9.5 

11.6 

- 2.4 

11.9 

33, 

294 

- 0.4 

3.4 

-17.3 

16.9 

32, 

289 

3.2 

4.0 

10.4 

-7.2 

32, 

290 

2.7 

4.0 

4.5 

-1.8 

32, 

291 

- 1.2 

- 0.2 

- 4.8 

3.6 

32, 

292 

18.7 

19.8 

10.5 

8.2 

32, 

293 

4.2 

5.7 

-6.8 

11.0 

32, 

294 

2.8 

2.6 

-13.5 

16.3 

31, 

289 

-25.1 

-24.2 

-20.8 

-4.3 

31, 

290 

-37.3 

-35.6 

-36.7 

-0.6 

31, 

291 

-26.1 

-24.3 

-31.1 

5.0 

31, 

292 

-17. 3 

-15.5 

-26.9 

9.6 

31, 

293 

- 7.1 

- 6.3 

-22.0 

14.9 

31, 

294 

- 5.7 

- 6.6 

-22.7 

17.0 

RMS 

11.5 

11.4 

13.8 

8.3 

RMS from first 
21 blocks 4. 8 

5.9 

6.7 

6.8 





















6.2 The Philippines Area 


Data : Adjusted altimeter data in the area 0° <(p < 40° and 120° < X < 160°. (For 
details see Rapp, 1977b.) 

In this area 16 5° mean anomalies and 5 l°x 1° anomalies were 
predicted. For the prediction of the 5° anomalies approximately 11300 observa- 
tions in 48 arc segments were included and the total computation time (IBM 370/168) 
was 1. 5 minutes. The results are shown in Tables 6. 5 and 6. 6. The RMS differences 
between predicted (weighted) and terrestrial (see Rapp, 1977a) anomalies was 7.7 
mgal. The corresponding difference between the GEM 9 and terrestrial anomalies 
was 11.2 mgal. 

The 5 l°x 1° blocks were particularly selected across a trench with 
very large variations of the terrestrial anomalies. . The results of these predictions 
are given in Table 6. 7. Large discrepancies are obtained between predicted and 
terrestrial anomalies. The main reason for the bad results is probably the 
following. Theoretically, the method requires that the profiles of observations 
are extended to infinity [formula (2.7)]. In practice, we can usually truncate the 
arcs at, say, 5° outside the prediction area without significant loss of information. 
However, in the present case the land masses (islands) on the west side of the 
trench truncate the arcs at a few degrees, which has’ a serious impact on the 
prediction results in such a disturbed area. 

For comparison, we report some l°x 1° anomaly predictions across 
the Mariana Trench, where there are no. disturbing land masses (Table 6.8). An 
obvious improvement of the GEM 9 anomalies is obtained, and the result is rather 
good if we consider the limited number of arcs. The maximum allowed distance 
arc-center of block seems to be very important in these very disturbed areas. 

See section 6,3, 


6.3 Dependence on the Number of Arcs 

In the previous computations, the maximum perpendicular distance 
from the center of a block to an altimeter profile was 3°. It is of interest to find 
out whether a change of this distance of truncation (not to be mixed with the length 
of the altimeter profile outside the prediction area) will change the prediction 
results significantly. The RMS differences between predicted mean anomalies 
and terrestrial anomalies for various such distances are given in Table 6.9. In 
general, the di ff erences are decreasing with increasing truncation distance, at 
least to 3° 5. However, in the Mariana Trench the result is opposite. We conclude 
that a smaller distance of truncation should be adopted in regions with large 
disturbances than in gentle areas. 
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Table 6.5 


5° Equal Area, Predictions in the Philippines Area 
Block Numbers and Terrestrial Anomalies (Ag T ± nu) Refer to Rapp (1977a) 

Units: mgal 


Block 

No. 

Coordinates 




Spread 

#of 

Arcs 

(n) 

AgT ± m T 

GEM 

9 

Wtd. 

Direct 

438 

30°-25°, 129°-135° 

-0.6 

-0.1 

3.4 

3.2 

6 

-8.4 ±2.8 

11.9 

439 

30°-25°, 135°-141° 

20.0 

19.9 

4.4 

4.6 

7 

20.9 ±2.8 

13.4 

440 

30°-25°, 141°-146° 

3.1 

6. 5 

4.5 

8.0 

7 

20.1 ±2.5 

10.5 

441 

30°-25°, 146°-152° 

-0.3 

1.2 

5.9 

2.5 

3 

0.6 ±2.8 

5.9 

503 

25°-20°, 129°-134° 

6.5 

5.9 

2.0 

1.9 

10 

7.8 ±2.7 

2.5 

504 

25°-20°, 134°-140° 

14.2 

14.8 

2.5 

3.5 

17 

10.6 ±2.7 

10.4 

505 

25°-20°, 140° -145° 

19.3 

18.3 

3.1 

10.5 

11 

20.1 ±2.9 

13.7 

506 

25°-20°, 145°-150° 

3.4 

4.1 

3.8 

8.5 

10 

-0.7 ±2.9 

9.6 

571 

20°-15°, 130°-136° 

5.3 

5.4 

1.6 

1.6 

17 

1.6 ±2.8 

4.9 

572 

20°-15°, 136°-141° 

13.4 

13.4 

2.6 

2.5 

14 

10.0 ±2.9 

7.8 

573 

20°-15°, 141°-146° 

20.0 

20.4 

3.5 

5.3 

12 

35.2 ±3.0 

9.5 

574 

20°-15°, 146°-151° 

6.2 

1.7 

3.1 

13.0 

11 

-5.3 ±3.5 

5.2 

640 

15°-10°, 129°-134° 

11.7 

11.7 


2.1 

16 

12.4 ±3.1 

15.3 

641 

15°-10°, 134°-139° 

12.0 

12.1 


2.9 

15 

10.3 ±3.2 

10.6 

642 

15°-10°, 139°~144° 

■ 4.8 

4.5 

3.0 

11.5 

13 

-8.7 ±3.3 

7.7 

643 

15°-10°, 144°-149° 

-8.2 

-12.9 

|g£ 

19.7 

11 

-8.3 ±2.8 

5.6 


* estimated prediction error from altimeter data 
+ = [E(Agi- Ag) 2 /np, where Ag t = prediction from arc i 
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Table 6. 6 


Differences Between Predicted, Terrestrial and GEM 9 Anomalies 

in the Philippines Area 
Units: mgal 


Block 

No. 

Weighted 

Terr. 

Direct 
Terr . 

— 

GEM 9 
Terr. 

Weighted 
GEM 9 

438 

7.8 

8.3 

20.3 

-12.5 

439 

- 0.9 

- 1.0 

- 7.5 

6. 6 

440 

-17.0 

-13.6 

- 9.6 

-7.5 

441 

- 0.9 

0.6 

5.3 

- 6.2 

503 

- 1.3 

- 1.9 

- 5.3 

4.1 

504 

3.6 

4.2 

- 0.2 

3.7 

505 

- 0.8 

- 1.8 

- 6.4 

5.6 

506 

4.1 

4.8 

10.3 

- 6.3 

571 

3.7 

3.8 

3.3 

0.4 

572 

3.4 

3.4 

-2.2 

5. 6 

573 

-15.2 

-14.8 

-25. 7 

10.5 

574 

11.5 

7.0 

10.5 

1.0 

640 

- 0.7 

- 0.7 

2.9 

- 3.6 

641 

1.7 

1.8 

0.3 

1.4 

642 

13.5 

13.2 

16.4 

- 2.9 

643 

0.1 

- 4.6 

13.9 

-13. 8 

RMS: 

7.7 

7.1 

11.2 

6.9 
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Table 6. 7 


Prediction of l°x 1° Anomalies in the Philippine Trench. 
Terrestrial Anomalies (Agr) According to Rapp (1977a). Units: mgal. 
Maximum Allowed Distance Arc-Center of Block: 3°. 


<p X 

NW 
Comer 

Pred. 

Means 

m A * 

4 * 

Spread 1 

# of 
Arcs 
<*) 

AgT 

±mr 

GEM 

9 

Wtd. 

Direct 


124° 

22 

19 

13 

8 

5 

70 ± 12' 

33 


125° 

31 

21 

8 

38 

9 

153 ± 20 

31 


126° 

15 

8 


28 

j mm 

-107 ± 25 

29 


127° 

20 

14 


22 


31 ±10 

27 


128° 

21 

18 


14 

wm 

43 ±10 

25 


* estimated prediction error from altimeter data 
t = [ S(Agj - Ag) 2 / n]& , where Ag t is the prediction from arc i. 


Table 6.8 

Prediction of l°x 1° Anomalies in the Mariana Trench. 
Terrestrial Anomalies According to Watts (.1975). Units: mgal. 
Maximum Distance Arc -Center of Block: 2°. 


<p X 

NW 
Corner 

Pred. 

Means 

m A * 

Spread^ 

#of 

Ares 

( n ) 

AgT 

±mt 

GEM 

9 

Wtd. 

Direct 

12°, 

139° 

6 

5 

16 

10 

3 


9 

12° , 

140° 

10 

10 

30 

- 

1 


8 

12°, 

141° 

- 12 

- 13 

23 

11 

2 

-172 ±25 

8 

12°, 

142° 

-111 

-117 

25 

44 

2 

-172 ±25 


12°, 

143° 

- 65 

- 61 

18 

75 

4 

-142 ±25 


11°, 

139° 

15 

15 

21 

0 

2 

18 ±15 


11°, 

140° 

- 25 

- 25 

31 

- 

1 

50 ±15 

9 

11°, 

141° 

2 

2 

23 

4 

2 

- 72 ±20 

8 

11°, 

142° 

4 

4 

24 

6 

2 

- 2 ±15 

MM 

11°, 

143° 

1 

3 

20 

19 

3 

15 ±10 

Kfl 


* estimated prediction error from altimeter data 

"i* = [ S (Ag j - Ag) 2 / n] s , where Ag i is the prediction of arc i , 
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Table 6.9 


RMS Differences (Weighted) predicted -Terrestrial Mean Anomalies 
for Various Distances of Truncation, 
i.e. Maximum Distances Arc - Center of Block 


Area 

Blocksize 

No. 

of 

Blocks 

Truncation Distance 

m 

m 

m 

m 

B 


4° 

Calibration 
1° x 1° 

30 


13.3 

12.2 

11.7 

11.5 

11.5 

11.7 

Calibration 
l°x 1° 

21 


■ 7.8 

6.2 

5.2 

.4.8 

H 

4.6 

Calibration 
’ 5° 

6 



4.1 

4.0 

3.8 

3.8 

4.1 

Philippines 

5° 

16 

10.0 

7.6 


n 

M 


■ 

Mariana Trench 
l°x 1° 

10 


69 

68 

75 

80 

81 

87 


Units : mgal 
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Finally, we compare 21 l°x 1° anomaly predictions for different 
numbers of included arcs (estimates). These numbers are selected simply by 
including only every 2-5 estimate of the original set of estimates. The results 
are given in Table 6. 10. 


Table 6. 10 

The Relation Between the Mean Number of Estimates (arcs)/ Block 
and the RMS Difference Prediction- Terrestrial Anomaly ( 1° x 1°) . 

Units: mgal. 


Mean No. of A res /block 

29.3 

14.9 

9.8 

7.4 

5.8 

RMS Difference 

4.6 

4.9 

4.8 

5.5 

6.7 


The table shows a trend of decreasing differences with increasing number of arcs. 
We can therefore expect that the inclusion of all available arcs in the area will 
improve the prediction results of sections 6. 1 and 6. 2. 


6.4 Test of Accuracy Estimates 

The tables of sections 6. 1 and 6.2 include the error estimates m A and mi 
of the estimates from altimeter and terrestrial data, respectively. As these 
different set of data are uneorrelated, we have 

cttot = E i { ( Ag A , - Ag T f } = E { mf } + E { m f } 
or 

(6.3) Or 0 T a = <7 A 3 + err 2 


where a 2 and cr 2 are the theoretical variances of altimeter and terrestrial data. 
Now we assume that the error of each Ag A and Agr estimate is normally 
distributed and that the variances are approximately constant in each prediction 
area. Then (6. 3) may be tested for the RMS estimates of Otot, of and of of each 
area, because (Ag A — AgT) 2 /(m A 2 + m 2 ) is approximately F - distributed with 
f - degrees of freedom (f being the number of blocks). This test is performed in 
Table 6.11, where F is the ratio of (Ag A - Ag T ) 3 and m A s + m 2 , with the largest 
outcome in the numerator. This ratio is compared with Ff t f fa /a> where a is the 
risk level. Formula (6. 3) is accepted if F< F f ,t,ct/a • The test shows 
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that the hypothesis (6. 3) is accepted in the calibration area, but is rejected in 
the Philippines area. Although the assumptions of constant prediction error within 
each area and uncorrelated estimates of various blocks are approximations , the 
test gives some idea of the validity of the estimated prediction errors. We 
conclude that only in the Philippines area (6. 3) is rejected, probably due to too 
small error estimates m A . 


Table 6.11 

Test of Formula (6. 3) 
Risk level (a.) = 5% 


Area 

#of 


F 

f v>0! / 2 

3m 

II 

Blocks ize 

blocks 

(f> 

Alt. 

Terr. 

m A 

m T 


Oa + of ? 

Calibration 
l°x 1° 

30 

11.5 

3.2 

13.6 

1.5 

2.1 

Yes 

Calibration 

5° 

6 

3.8 

1.0 

2.4 

2.1 

5.8 

Yes 

Philippine 

5° 

16 

7.7 

3.3 

2.9 

3.1 

2.8 

No 


7. Conclusions and Final Remarks 

In this report we have shown that if the data are given along profiles 
(ares) an approximate version of least squares collocation may be performed in 
the frequency domain rather than in the direct domain. In this way the auto- 
covariance matrix, that has to be inverted, is transformed to diagonal form. In 
the present application, each profile is processed separately and the estimates 
from all arcs are finally joined in mean estimates. A considerable gain in 
computer time is obtained in comparison with conventional collocation technique. 
According to Thomas (1977) the necessary computer time for n observations is 
reduced from ~ n 3 to ~n log n by using the Fourier transform method. For 
example, approximately 11300 observations in 48 arc segments were included • 
in one of the'present computations, 16 5° anomalies were predicted in 1, 5 
minutes computer time (IBM 370/168). In addition, new information that becomes 
available can easily be included in the predictions. 
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One should also consider that numerical difficulties occur in the 
direct method when the positions of two or more observations are close. Thus 
a filtering of the data is necessary, which probably implies a loss of information. 
The Fourier transform technique does not suffer from such limitations. 

A disadvantage of the method is the requirement of equal spaced 
data along each arc. As the altimeter observations do not have a regular spacing 
along the profile (but do in time), such an approximation might be. erroneous for 
point estimation, but will hardly effect the mean block predictions. 

A summary of the computations is given in Table 7. 1. The altimeter 
predictions are weighted means. 


Table 7. 1 

Summary of Prediction Results 
Approximate Numbers of Used Arcs and Observations 
Computer: IMB 370/168 


Area 

No. 

No. 

No. 

m 

RMS Diff. 

m A * 

^^9 

CPU 

Blocks 

of 

of 

of 


(mj 

gal) 

(mgal) 

Time 


Obs. 

Arcs 

Blks 


GEM 

Terr. 

Pred. 

Terr. 



(min.) 

Philippine 

5° 

11300 

48 

16 

60 

10.5 

7.7 

3.3 

2.9 

1.5 

Calibration 

5° 

12800 

85 

6 

60 

7.1 

3.7 

1.0 

2.4 

2.1 

Calibration 

1° 

6550 

41 

30 

75 

13.8 

11.5 

3.2 

13.6 

8.0 

Calibration 

1° 

6550 

41 

21 

75 

6.7 

■ 4.8 

3.1 

12.4 



* estimated prediction error from altimeter data 
"f" estimated prediction error from terrestrial data 
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The table shows that a reasonable improvement of the GEM 9 anomalies 
is obtained by adding the altimeter information. However, this improvement is less 
obvious for the 30 l°x 1° anomalies, because of the large errors of the terrestrial 
data. In section 6. 3 we have demonstrated that the inclusion of more arcs in the 
predictions will probably give even better results. Due to the simplicity of the 
method, we can easily include all available arcs in the area. As in the direct 
collocation method, mean geo id undulations may be estimated by substituting the 
covariance function Cakn by Cnn* 

We feel that the estimated prediction errors are somewhat too optimistic. 
However, only for the 6 5° equal area blocks in the Philippines area we found some 
evidence that this is the case. 

Finally, it should be stated that the method fails in coast (or island) 
areas with large disturbances. The reason is probably the short truncation distances 
of the arcs in the direction of the land masses. Further investigations are needed to 
reveal the impact on the final predictions of the truncation distances (both with respect 
to the length of the arc and with respect to the distance arc - predicted block), the 
number of discrete intervals (m) in the numerical integrations and the choice of 
covariance functions. 
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Appendix 


We are going to show that the weight: 


(A.l) 


ho = ir (h© + hi) 


according to Moritz (1967, p. 49) is incorrect and that the correct formula is: 


(A.l) 


ho - l/4 (ho + hq ) 


For the proof we need the following identity given by Blackman and Tukey (1958, p, 80) 


^ cos (0 + 2hu) = 


sin {(a+b+l)u} 


sm u 


cos {0+ (b-a)u] ; u/0 


Inserting a — -1, b-m-1, 0=0 and u = 4rr/2m, we arrive at the following formula: 
(A. 2) 


Furthermore 
(A. 3) 


V „ „ , sin {(m-l)£rr/2m} jCtt 

) cos (kArr/m) = - - J cos = 

^ sin (Air /2m) . 2 


k=l 


{ 


(- 1 ) 


l/z 


for 1= 1, 3, 5, ... 

for £= 2, 4, 6, ... 


a-i 

^ cos (kjgrr/m) = m - 1 for i= 0 


k = l 


The relations (A. 2) and (A. 3) will be used in the following derivation of the weight h 0 . 
From (3. 8) or Moritz (ibid.) we have: 

Aco 


a- 1 


✓ iico r v , -1 

h 1 - Ho + 2 2^ H q cos (qj£n/m) + H„ cos ( j&tt ) 


1=1 
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In the special case of pure prediction (of X from the true values of X) all system 
functions H a are unity. Thus we obtain: 

a— 1 

h/= c[l+2^ cos (q^TT/m) + cos (£tt) j 

<i=i 


where [from formula (3.3)} : 


c = Ato/2rr = l/2 m At 


Inserting (A. 2) and (A. 3) we arrive at; 


(A. 4) 



[l+(-l) 1/2 ]/m 


{, = 0 

& = If 3, 5} . • • 
£- 2, 4, 6, ... 


Next we determine the smoothed weights from the raw values {see formula (3.9) 
and Moritz (ibid.) } . However, we use the modification: 


(A. 5) 


ho = P (h 0 " + h/) 


where p is a constant to be determined in such a way that the prediction by formula 
(3. 1) is correct (without error), i. e. 


(A. 6) 


D + l 

X a =Aty hi (X k . } + X k+} ) = X k 

c i 

S— o 


For simplicity we assume that: 

X k = X = constant for all k 


Furthermore, we assume that m is so large that (A. 4) becomes (without loss of 
significance) : 


(A. 7) 



&= 0 


otherwise 
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Then the smoothed weights are given by (3.9) with the modification (A. 5).. 
Inserting (A. 7) we obtain: 

h 0 = p/At , h x =1/4 At 

and h q = 0 for q > 1 


Combining these weights with (A. 6), we obtain the following equation: 


X(P+P+ 1/4 + 1/4)= X 


with the solution for p: 


Thus we have proved that: 


p = l/4 


ho — 1/4 (ho + h^) 


is the correct expression for ho . 
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The Determination of Gravity Anomalies from Geo id Heights 
• Using Least Squares Collocation* 


by 


Richard H. Rapp 
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the title "Anomalies Recovered from Geos-3 Altimeter 
Data Using Least Squares Collocation" 
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Abstract 


Geos-3 altimeter data provided by the Wallops Flight Center has been 
processed using techniques developed by Rummel to determine consistent geoid 
height information without need for a precise Geos-3 orbit. Data was selected 
from fifty-three processed arcs to determine three 5° equal area anomalies, 
two 2°x 2° anomalies and two l°x 1° anomalies. These blocks were located 
so as to maximize the amount of altimeter data in and surrounding the block 
whose anomaly was to be estimated. The anomalies were estimated using the 
technique of least-squares collocation using both an ellipsoidal reference field 
and a reference field defined by the GEM 7 potential coefficients taken to degree 
12. The predicted anomalies showed good agreement with estimates based on ■ 
terrestrial gravity data. The average predicted standard deviations were ±3 
mgals for the 5° blocks; ± 5 mgals for the 2°x 2° blocks; and ±10 mgals for 
the l°x 1° blocks. 


- 52 - 



Introduction 


For several years discussions have taken place on the recovery of ' 
gravity anomalies from satellite altimeter data. These discussions have been 
limited to simulation studies or studies using data other than altimeter data. 

In a paper presented last year at the American Geophysical Union Meeting, I 
described a series of simulation studies' that showed the expected accuracy of 
the recovery of various mean gravity anomalies from altimeter data postulated 
to be of varying density and of various accuracy. At- this time we now have 
sufficient real altimeter data from Geos-3 to carry out actual anomaly recovery 
with the results being compared to actual estimates based on terrestrial gra- 
vity measurements. This paper reviews the theory, the data used, the results 
and gives recommendations for future work. 


The Theory 


The general method that we use is called least-squares collocation - 
(Moritz, 1972). The application of this general method to the use of altimeter 
data for anomaly recovery work has been discussed by Smith (1974) and Rapp 
(1974) among othe rs . 

We will use the simplified model of collocation where we assume that 
the observations do not contain systematic components that are described by a 
term of the form AX where X is a vector of parameters of unknowns. In this 
case the basic equation of least-squares -collocation is; 


( 1 ) 


x = s r + n 


where: 


x is a measurement (which in our case 
is a geoid undulation derived from the 
measured altimeter data); 
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is the random signal part of x; 
n is the random noise of the observation. 


One result of least squares collocation can be the estimation of a 
signal, s, of a quantity that we desire to determine, which in our case is 
gravity anomalies, Ag. The collocation solution of (1) for a predicted signal 
is: 


( 2 ) 


— C i »t (C*f *? + C nn ) £ 


where; 


C*»* is the covariance between the signal to 
be estimated and the signal component of 
the measurement; 

C,v is the covariance between the signal com- 
ponents of the measured data; 

C na is the error variance-covariance of the 
measurements. 


The variance-covariance matrix of the predicted signal is: 


( 3 ) 


Ea» — C a « - Cg»’(C*V + C n a) ^ Cb ! » 


where C«* is the covariance matrix for the si gnal s being estimated. 

For the prediction of a single gravity anomaly at point P from 
a set of altimeter data, a, we can write equations (2) and (3) in the fol- 
lowing form: 
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( 4 ) 


Agp = (C aa + C^- 1 a 


(®) m P - C Sp gp - Cg^d {Caa + C nn ) 

where rrjj, is the standard deviation of the prediction and C K t is the expected 
mean square value of the anomaly block being estimated. 


The Covariances 


hi order to' evaluate equations (4) and (5) it is necessary to compute 
the covariances between the altimeter (undulation) data, as well as the cross- 
covariance between the anomaly and undulation data. One might consider the 
use of covariances derived from actual data in the area of computation. The 
determination of such covariances requires sufficient data which in most cases 
is lacking. We consequently use the concept of a global covariance function 
which is isotropic and stationary and which can be computed analytically once 
certain parameters of an anomaly degree variance model are specified. The 
details of the formulation of such models and their derivations .leading to a 
computer subroutine that can be used for the generation of the needed covari- 
ances are given in Tscheraing and Rapp (1974). All covariances used in this 
paper were computed with subroutine COVA given in that report. 


The covariances that are given by COVA are regarded as covariances 
between point quantities. However, we are estimating mean anomalies, not 
point anomalies so that special consideration must be given to derive the appro- 
priate covariance value for final use in equations (4) and (5). Thus, for ex- 
ample, with Ag p considered a mean anomaly, Cg pa would be the covariance be- 
tween the mean anomaly block and the altimeter measurement (generally re- 
garded as a point measurement) . These mean covariances or variances can be 
computed by the numerical integration of the point covariance functions over the 
block in question (Heiskanen and Moritz, 1967, p.277). 


- 55 - 



The Reference Field 


Anomalies or other gravimetric dependent quantities can be given with 
respect to an ellipsoid reference field or to a higher degree field defined by a set 
of potential coefficients given to some specified degree, N. It might be ejected 
that more accurate predictions could be carried out by using a high degree refer- 
ence field but this expectation must be judged against the accuracy of that refer- 
ence field. 

In this paper, for comparison purposes, we will carry out predictions 
based on the usual ellipsoidal model, and then with a reference model defined by 
the GEM 7 potential coefficients (Lerch et als. , 1975) taken to degree 12. hi all 
cases the final predicted anomaly is to be given with respect to an ellipsoidal re- 
ference model. 

To carry out this prediction using the higher degree field we first sub- 
stract from our altimeter data (actually the undulation implied by the altimeter 
data), the undulation implied by the potential coefficients. We then predict the 
mean anomaly with respect to the reference field using the covariances that have 
had the reference field components removed. To obtain the final predicted ano- 
maly we then add the anomaly in that block being predicted implied by the poten- 
tial coefficients .adopted for the reference field. 

In spherical approximation, equations used for computing the refer- 
ence undulations and anomalies are; 




( 6 ) Agj, = G ^ (n - 1) ^ (C^eosmX + S nB sinmX) P ftn (sin <0 

n=£ m=0 


( 7 ) 


N n 


N H = R ^ ^ (C^eos mX + S 3 a sinmX) P Ra (diri^). 


n =3 b = o 


where; 



C na , S ns . . are fully normalized potential coefficients; 

3.0 ~ C2,0 ** C2,c(rs?); 

C 4 *o = C 4 ,o ~ C 4>3 ^ r9 ^; 

P DB ...... are the folly normalized associated Legendre functions 

<£>, X geocentric latitude and longitude 

R, G average values of the earth's radius and gravity. 


The Data and Its Accuracy 


The Geos-3 data being analyzed in this study represents data received 
from the Wallops Flight Center for the time period April 21 through May 20, 1975. 
Of the data analyzed in this paper 48 arcs were received in February 1976, while 
5 arcs were received in October 1975. This data consisted both of the long pulse 
and short pulse mode with the angular spacing between data points being 0? 12 for 
the long pulse mode and 0?19 for the short pulse mode. 

The actual altimeter data, converted to raw undulation data, is contam- 
inated by several errors such as orbit error, altimeter bias etc. It is thus neces- 
sary to process the data that we have received to remove such errors and to put 
the data in a consistent system, .This processing is described in a paper by Rum- 
mel (1976). The input data to the coHocation prediction are the unfiltered geoid 
undulations obtained by Hummel at each of the altimeter data points. These undu- 
lations however, may still have errors (other than observation errors) in them 
which arise from sea state effects and unmodeled errors in the original proces- 
sing at Wallops or in the error model adopted.by Rummel. These errors are ex- 
pected to be small. 

In order to evaluate the C nn matrix in equation. (4) and (5) it is necessary 
to have an accuracy estimate for each measurement. Such an estimate was avail- 
able from the data supplied by Wallops. Thus estimates varied from approxi- 
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mately ±0. 5m to ±8. 5 m with average values being on the order of 1. 2 m 
(Rummel, 1976). Preliminary anomaly estimations were done assuming 
uniform accuracies of ±lm, an error estimate for each track, and the indi- 
vidual error estimates. The predicted anomalies showed only a slight de- 
pendence on the accuracy assumed for the data. Consequently, we assumed 
that the accuracy of the undulations derived from the altimeter data were 
those given for each altimeter measurement through the processing proce-. 
dure at the Wallops Flight Center. 


Anomaly Predictions 


We chose to estimate three 5° equal area anomalies, two 2°x 2° 
anomalies, and two l°x 1° anomalies. The location of these anomalies was 
chosen to maximize the altimeter tracks in and around the blocks. The lo- 
cation of these blocks and an alphabetic block designator are shown in Figure 
One, with the block coordinates being given in Table One. 


Table One. Coordinates of Blocks Used in the Anomaly 
Prediction. 


Letter 

Designation 

Block 

Size 

<p 

X 

A 

5°E.A. 

45° - 40° 

299°- 306° 

B 

5°E.A. 

40°- 35° 

291°- 297° 

C 

5°E. A. 

35° - 30° 

295°- 301° 

D 

' 2° x 2° 

40° - 38° 

291°- 293° 

E 

2° x 2° 

34° - 32° 

296°- 298° 

F 

l°x 1° 

39° - 38° 

291°- 292° 

G 

1° X 1° 

34° - 33° 

296°- 297° 
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The gravity data used to obtain the check anomalies was an updated 
version of the DMA AC (1973) l°x 1° data. This updated version is described 
in Rapp (1975). The l°x 1° data is taken directly from the tape. The 2°x 2° 
and 5° equal area anomalies were formed by a direct mean of the l°x 1° ano- 
malies in the block. The anomalies and their standard deviations are given in 
Table Two. 

In carrying out the actual predictions several decisions have to be 
made. Some of these areas of concern are: 

1 . Reference Model 


We have previously discussed the use of an ellipsoidal reference 
model and a model defined by the GEM 7 coefficients taken to degree 12. 
Computations have been made with both fields to compare results. Other re- 
ference fields and different maximum degrees could have been tested but were 
not for this paper. 

2. Covariances 


The covariances used were the global covariances obtained from 
COVA as previously explained. No optimization of the covariances was done 
to fit the data in the test area. This may be considered in the future. 

3. Data Cap Size 

Data within a cap about the center of the block being predicted is 
to be chosen. The cap is chosen large enough so that data outside this cap does 
not influence the predicted anomaly or its accuracy within some specific toler- 
ance level. Simulation studies indicated cap sizes on the order of 5° to 7° 
were reasonable when dealing with an ellipsoidal reference model while a cap 
size on the order of 3° to 4° was reasonable when a degree 12 field was used. 
In this paper results for both cap sizes will be shown for the 2°x 2° and l°x 1° 
predictions. 

4. Data Density 

Although the use of all known data points within the chosen cap 
might be considered for use in the anomaly prediction there are two reasons for 
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not using all data. First, the size of the matrix to be inverted in (4) and (5) 
is equal to the number of observations. Thus a large number of observations 
will lead to large matrices with potentially large inversion times. 

Second, the use of many points, close together, may yield some 
instability in the inversion as the elements of will become similar. 

We thus must choose our observations with some care. In com- 
putations made for this paper we selected only a certain number of. points along 
a track that fell within the specified cap. We did tests choosing every 2nd 
point, 3 rd point and in some cases 4th point. In some cases the choice of every 
other point led to some instability relative to the case of using every 3rd point. 
For the predictions presented in this paper every third point along a track was 
selected. Results for other point selections are available but they are not pre- 
sented here* 


Results 


The anomaly predictions are shown in Table Two. In the case of 
the use of the reference field to degree 12, the standard deviation of the predicted 
anomaly has been given assuming the reference field is known perfectly. 

Table Two. Anomaly Prediction Results. 


Block 

Size 

Terrestrial 
Anomaly (mgals) 

Predicted 

E* 

Anomaly 

N=12t 

Num. of 
Points 

A 

5°E. A. 

-17 ±3 

-16 ±6 

-20 ±3 

198 

B 

5°E.A. 

-24 ±2 

-24 ±5 

-25 ±3 

225 

C 

5°E. A. 

- 8±3 

- 3 ±6 

-11 ±3 

164 

D 

2° x 2° 

-31 ±7 

-33 ±7 

-32 ±5 

161 

E 

2° x 2° 

- 8 ± 7 

- 2 ±7 

- 8 ±5 

142 

F 

l°xl° 

-25 ± 14 

-30 ±11 

-29 £9 

165 

G 

l°x 1° 

- 8 ±15 

- 3 ±10 

- 9 ±8 

145 


* Prediction made with the ellipsoidal reference model, 
t Prediction made with GEM 7 to degree 12 as the reference model. 

9/1/76 
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With the exception of one or two cases the results using the ellipsoidal 
and GEM 7 reference models agree quite well, in addition the predicted anoma- 
lies agree well with the terrestrial estimates. For those tests using two different 
cap sizes the larger cap generally gave the better result. 


Conclusions 


The results obtained in this paper indicate that terrestrial anomalies 
can be recovered from the Geos-3 altimeter data using least squares collocations 
The predicted values agree with the known values within the range of the predicted 
standard deviations. These standard deviations for 5° anomalies was about ± 6 
mgal which were on the order of ±3 mgal in ideal circumstances. Nevertheless, 
a good percentage of existing oceanic 5° equal area anomalies are not known to 
±7 mgal so that the altimeter data can provide a real improvement to their deter- 
mination. 


The predicted accuracy of the l°x 1° anomalies was on the order of 
± 10 mgals (although the actual predictions were better than the standard devia- 
tion implied). Denser data should yield accuracies on the order of ± 5 mgals. 

Finally we note that these results indicate that the processing of the 
semi- raw altimeter data yields undulation information of high quality without 
the need for precise satellite orbits. 

Future work on anomaly determination awaits the release of more 
Geos-3 data. 
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